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Evolution of highly buoyant thermals in a stratified layer 
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Abstract. The buoyant rise of thermals (i.e. bubbles of enhanced entropy, but initially in pressure equilibrium) is 
investigated numerically in three dimensions for the case of an adiabatically stratified layer covering 6-9 pressure 
scale heights. It is found that these bubbles can travel to large heights before being braked by the excess pressure 
that builds up in order to drive the gas sideways in the head of the bubble. Until this happens, the momentum of 
the bubble grows as described by the time-integrated buoyancy force. This validates the simple theory of bubble 
dynamics whereby the mass entrainment of the bubble provides an effective braking force well before the bubble 
stops ascending. This is quantified by an entrainment parameter alpha which is calculated from the simulations 
and is found to be in good agreement with the experimental measurements. This work is discussed in the context 
of contact binaries whose secondaries could be subject to dissipative heating in the outermost layers. 
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1. Introduction 

Highly buoyant bubbles with large specific entropy ex- 
cess relative to the surroundings have been invoked by 
Hazlehurst (1985) in an attempt to explain the almost 
equal effective temperatures of the two components of con- 
tact binaries. 

As noted by Sinjab et al. (1990) there are strong paral- 
lels between the highly buoyant bubbles of Hazlehurst and 
the 'thermals' found to occur in the earth's atmosphere. 
Subsequently, Hazlehurst (1990) confirmed the existence 
of a formal relationship between Turner's (1963) treat- 
ment of thermals, involving entrainment of matter, and his 
own treatment of highly buoyant bubbles ('interlopers') in 
contact binaries, as these bubbles annex new material. 

In this paper we shall discuss the question of the va- 
lidity, from a fluid-dynamical standpoint, of the simple 
bubble or thermal picture. Wc shall then go on to show 
how it is possible to determine numerically the value of 
the entrainment coefficient a (here called av) which enters 
Turner's and several other investigations; the determina- 
tion of a related coefficient (here called am) entering the 
Hazlehurst theory is also discussed. 

We believe this to be the first attempt to evaluate the 
(previously semi-empirical) entrainment coefficient a on a 
fluid-dynamical basis. 



2. Model setup 

We adopt a basic setup of our model that is similar 
to that used normally to study convection in a strati- 
fied plane-parallel layer between impenetrable boundaries 
(e.g. Hurlburt et al. 1984). In particular, we use stress- 
free boundary conditions at the top and bottom, with a 
prescribed flux F at the bottom and a prescribed tem- 
perature Ttop at the top. Here, however, we assume the 
thermal equilibrium stratification to be adiabatic, so it is 
marginally stable to the onset of convection. Thus, when 
we insert a hot (buoyant) bubble, it will rise unaffected by 
the stratification (except for effects related to the growth 
and expansion of the ascending bubble), so there is no 
restoring force acting on the bubble as it rises. 



2.1. Adiabatic stratification 

Hydrostatic equilibrium requires that the weight of the 
atmosphere is balanced by the pressure gradient. However, 
if the entropy is constant, the pressure gradient can be 
written as Vp = V/i, where p is pressure, p is density, 
and h is enthalpy. We adopt a perfect gas for which h = 
CpT, where Cp is the specific heat at constant pressure. 
For constant gravitational acceleration .g > 0, this implies 
a linear temperature stratification that is given by 



^(cpr) = g, 



(1) 
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where z is depth, which is assumed to increase downwards. 
Thus, the vertical temperature profile of the basic state is 
given by 



T = Ttop + {z - ztop)g/cp 



(2) 



In the absence of any motion there is only radiative flux, 
F, for which we adopt the diffusion approximation, so 



-KVT, 



(3) 



where K is the radiative conductivity. Thermal equilib- 
rium requires V ■ = 0, so the 2;-component of the flux 
is constant but, because the temperature gradient is con- 
stant, this is only possible ii K = const. We now need to 
discuss the choice of K and other parameters. 

2.2. Choice of parameters 

We adopt nondimcnsional units by defining a unit length 
d, and our bubble will usually have the initial radius 
Ro = d (although we present initially some cases where 
i?o = 0.5 c?. We measure time in units of {d/gY^^, density 
in units of po (we choose p = po at the location where 
the centre of the bubble will be introduced), and specific 
entropy, s, in units of Cp. This corresponds to setting 



= 5 = Po = Cp = 1- 



(4) 



A relevant nondimcnsional quantity is the normalised in- 
put flux 



F 



(nondimcnsional input flux). 



(5) 



For secondaries of contact binaries this ratio is around 
10~^. Specifying F fixes K = CpF/g. There is however 
the numerical constraint that the mesh Peclet number, 
based on the sound speed Cg and the mesh width Ax, 



Pe 



'grid 



(6) 



should not exceed a certain empirical upper limit of 10- 
100. This limit arises from the fact that the advection of 
sharp structures tends to generate small ripples on the 
scale of the mesh which need to be damped. Here, Cg — 
IP/p-: 7 = Cp/cv is the ratio of specific heats, and x = 
K/ pcp is the radiative diffusivity. In the present paper we 
shall consider values of F in the range 0.001-0.005. In all 
cases we assume 7 = 5/3. 

We adopt cartesian coordinates (x, y, z) where z points 
downwards. The bubble centre is placed initially at a; = 
J/ = z = 0. In order that the bubble be initially sufficiently 
far away from the bottom boundary and that it can rise 
over a distance of at least a few times its own radius (which 
is unity in most of the models) , we chose the vertical extent 
of the computational box to be from ztop = —6 to Zbot = 2. 

Next, we fix the temperature stratification within the 
box by specifying the value of Ttop, which is held constant 
by the boundary condition chosen. The choice of this pa- 
rameter is restricted by numerical considerations. Since 
temperature is proportional to the pressure scale height. 



small values of Ttop iniply a short pressure scale height 
at the surface. However, numerical stability and accuracy 
considerations require that the scale height cannot be less 
than at least 2-3 mesh zones. The non-dimensional pres- 
sure scale height at the top is 



Co = Hp^top/d = (cp - Cv)Ttop/gd, 



(7) 



where we have included d and g factors, even though they 
are unity. For runs with Nz = 50 meshpoints in the ver- 
tical direction, we were able to use — 0-3. For orienta- 
tion, we note that this yields a ln(pbot/?'top) of 6.2 pres- 
sure scale heights between top and bottom of the box. 
The local pressure scale height at z = is then given by 
0.4 X jztopl + Co = 2.7. Using Eqs (||) and (0) we also see 
that the top of the adiabatic atmosphere would be at 



ztop - 2.5^0 (for 7 = 5/3), 



(8) 



but this should be outside the computational box. For 
ztop = —6 and Co = 0.3 this gives Zqo = —6.75. 

For a perfect gas we have p/p = (cp — Cv)r and 
s = Cvlnp — Cplnp (in dimensional form). The density 
stratification follows then from hydrostatic balance and 
the assumption that the entropy of the imperturbed model 
is constant: 



p/po = {i-z/zo,y/<^^~'i 



(9) 



We mentioned already that we express density in units 
of Po, which is the density at z = 0, which is where the 
bubble will be positioned. Using this together with the 
definition of s we can use Eq. (H) to express the value of 
the background entropy. 



So = 0.6 ln(-0.4zoo) (for 7 = 5/3). 



(10) 



For Zoo = —6.75 this gives sq — 0.596. 

This completes the definition of the background state. 
We now turn to the discussion of the model equations and 
initial and boundary conditions adopted. 

2.3. Governing equations 

In the dynamical case the specific entropy is not only af- 
fected by the radiative flux divergence, but also by the 
rate of viscous dissipation, so 



pT^ = V • KVT + 2i^pS^, 



(11) 



where D/Di ~ d/dt + m • V is the lagrangian derivative, 
1/ = const is the kinematic viscosity and 

Sjj = ^{djUi + diu-j - ^SijdkUk) (12) 



is the (traceless) rate of strain tensor. Equation ( |ll| ) is 
solved together with the momentum equation, 

1, 



Du 1„ 

ut p 



g+-V ■{2vpS), 
P 



and the continuity equation 
D Inp 



= -V • u. 



(13) 



(14) 
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We solve Eqs @ and using the sixth order 
compact derivative scheme of Lele (1992) and a third order 
Hyman scheme for the time step. For earlier applications 
of this code see Nordlund & Stein (1990) and Brandenburg 
et al. (1996). 

The value of v is dictated again by numerical consid- 
erations, and in practice we take v ~ but note that ly 
is independent of z whilst x is not. In all cases considered 
below we have used — 6 x 10^^. 

2.4. Initial and boundary conditions 

The bubble centre is placed initially atx = y = z = and 
has an entropy profile of the form 

s = so + As$(r), (15) 
where 

$(0 = I if^<^" (16) 

[ otherwise 

is a profile function defining the initial shape of the bub- 
ble. Here r is the initial distance from the centre. Initial 
pressure equilibrium requires that the increase of s is com- 
pensated by a corresponding decrease of In p, so 

lnp= ln( 1- — ) - As$(r). (17) 

7-1 V ^oo / 

The initial entropy excess of the bubble is given by the 
parameter As. [In all cases presented we take As — 0.5, 
which corresponds to the value used by Hazlehurst (1985), 
who adopts units where the nondimcnsional specific en- 
tropy is larger by a factor of 5. We note that larger values 
of As make the bubble rise faster, but we found that even 
for As = 2 the motion remained subsonic] 

We assume stress-free boundaries at the top and bot- 
tom, and prescribe the value of T at the top and dT/dz 
at the bottom, as is usual for convection calculations (e.g. 
Hurlburt et al. 1984). We use periodic boundary condi- 
tions in the x and y directions. The horizontal extent of the 
box, I a: I < and \y\ < Ly, is varied between — Ly ~ 4 
and 16. 

2.5. Allowing for three-dimensional effects 

In order to assess the fragility of the bubble during its 
ascent we have adopted in many cases substantial ini- 
tial velocity perturbations. In Fig. |l| we show a three- 
dimensional representation of the entropy for a run 
with As = 0.5 using initial velocity perturbations with 
Wmax/cs = 0.4 and Wims/cs = 0.02. As is evident from 
Fig. the entropy of the blob is hardly affected by these 
perturbations and only near the surface does one see 
strong perturbations. 

It is interesting that the bubble remains an entity dur- 
ing much of its ascent. In fact, even when the initial con- 
dition is quite different, bubble- like structures tend to de- 
velop. As an example we show in Fig. |^ a case where we 




Fig. 1. Three-dimensional representation of the specific 
entropy. The initial velocity perturbations are noticeable 
mostly near the top layers, but the bubble itself remains 
fairly axisymmetric. 100 x 100 x 50 meshpoints. 



have introduced an almost uniform horizontal layer of en- 
hanced specific entropy with a gaussian vertical profile 
initially. We have superimposed random small scale per- 
turbations to get the buoyancy instability started. 




Fig. 2. As Fig. |l|, but for the case of a horizontal layer 
of enhanced specific entropy initially. 50 x 50 x 100 mesh- 
points. 
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In the following we study in more detail the dynamics 
of an isolated buoyant bubble or thermal, as it shall also 
be referred to. 



3. Dynamics of isolated thermals 

In the following we consider cases with different degrees 
of stratification and different extents of the computational 
domain. 



t = 0.0 



t = 5.0 





3.1. Modest stratification 

We begin by considering first vertical cross-sections of en- 
tropy and velocity; see Fig. |^. The archimedian buoyancy 
force is largest in the middle of the bubble, and that is also 
where the vertical velocity is largest. On both sides of the 
bubble the velocity turns over, as expected (compare with 
observations of thermals described by Scorer 1957). 

A rather different impression is obtained when look- 
ing at the bubble in a comoving frame of reference; see 
Fig. ^. In this frame there is a stagnation point and hence 
there is a clear distinction of regions inside and outside the 
bubble. The flow pattern generally conforms with the no- 
tion of bubbles behaving like balloons with a more-or-less 
well defined surface and gas flowing around this surface. 
However, the bubble clearly grows in size and even its 
mass grows during its ascent. 

It should be noted that this addition of new material 
with 'different' specific entropy does not contradict the 
bubble concept. This is because there is no entropy dis- 
continuity between the bubble and the surroundings, so 
that the entropy of a captured particle can be changed 
gradually - by friction and by thermal diffusion - as it 
enters the bubble. 

By the time i = 30 the effects of the lateral boundaries 
have begun to affect the evolution of the bubble. Therefore 
we show in Fig. || the case of a wider box {L^ = 8). Note 
that there is now a noticeable flow speed even beyond 
I a; I =4 (the extent of the box in the previous case). In 
this calculation we have also included strong velocity per- 
turbations, but the overall flow pattern is still dominated 
by the rising bubble. 

In order to quantify the rise and the growth of the 
bubble in detail, we define the bubble B as all points in 
space where s > Scrit, with Scrit just a little larger than the 
background value, which is here sq — 0.596, so we chose 
Scrit — 0.597. The volume of the bubble is then estimated 
as 



t = 10.0 



t = 20.0 



V{t) ^ / dV, 

J B 

and the volume radius of the bubble is 

R{t)^{V/^^Y'\ 

The mass of the bubble is 



(18) 



(19) 




-3 



t = 30.0 




-a 




-2 3 4 



Fig. 3. Velocity vectors superimposed on a grey scale 
representation of the entropy (dark indicates high entropy; 
all panels have the same grey scale). The velocity is shown 
in a fixed frame of reference. The initial specific entropy 
excess of the bubble is As = 0.5 and its initial radius is 
-Rq = 0-5- The single contour shows the position where 
s = Scrit = 0.001 + sq. F ^ 0.005, 50^ meshpoints. 

In Fig. ^ we plot R(t) and M{t). Both functions increase 
monotonically, except for some minor departures at late 
times when the bubble has reached the top of the layer. 
Note also that at early times (t < 4) i? increases some- 
what faster than at later times. Qualitatively this type of 
behaviour is expected, because radiative diffusion causes 
structures to grow proportional to t^/^, which causes an 
infinite slope of i?(t) at t = 0. 

In order to make further comparison with Hazlehurst 's 
theory of buoyant bubbles, we measure the position of the 
centre of mass of the bubble, 



^bubble 



zpdV / I pdV. 



(21) 



M(t) 



pdV. 



Since z decreases upwards we define the height of the bub- 
ble as h = — Zbubbie- The height h and velocity v — Ah/ At 
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t = 0.0 



t = 5.0 




-2 



t = 10,0 




-3 



t = 20.0 





-5 



Fig. 5. Same as Fig. ||, but for i = 30 and for a wider 
box {Lx — 8). Although in this calculation strong three- 
dimensional perturbations have been applied initiaUy, the 
effects on the shape of the bubble are negligible. F = 
0.005, f 00 X f 00 X 50 meshpoints, As = 0.5, Rq = 0.5. 



t = 30.0 




t = 40.0 





10 



20 
t 



30 



40 



Fig. 4. Same as Fig. but now the velocity vectors are 
shown in a comoving frame of reference. 

of the bubble are shown in Fig. 0. Note that the height 
seems to approach a maximum near h — 3.5, so the centre 
of mass of the bubble does note quite reach the top of the 
box. (Below we shall show that for larger bubbles, i?o = 1, 
the maximum height is even less, suggesting that this is 
at least partly a geometrical effect; we shall also see that 
the top of the bubble does reach the top of the box in all 
cases.) 

The momentum of the bubble, Afw, is plotted in Fig. || 
and compared with the time integrated buoyancy force, 

f{Md-M)gdt, (22) 
Jo 

where 

Md{t) ^ [ pdV (23) 
Jb 

is the mass of displaced material and p is the density of the 
undisturbed medium. Between t = 5 and 15 the momen- 
tum of the bubble is somewhat larger than expected from 
the buoyancy force. This discrepancy depends somewhat 




Fig. 6. Radius and mass of the bubble shown in Fig. |[ 
F = 0.005, 50^ meshpoints. As = 0.5, Rq = 0.5. 

on the definition of the boundary of the bubble. However, 
more dramatic is the sudden loss of momentum of the bub- 
ble after t « 20, whilst the buoyancy force, as estimated 
by Eq. (p2|), continues to operate beyond this time. 

There are at least two possible reasons for this sud- 
den braking effect. One reason could be that the blob 
gets too close to the top and loses momentum simply be- 
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mechanism. Thus, we now turn to the first possibiUty, of 
which we can distinguish two variants. It is possible that 
the pressure build-up near the top is either an artifact of 
the top boundary being impenetrable, or it could be sim- 
ply a feature of strong density stratification which causes 
the bubble to expand rapidly sideways. In order to drive 
strong sideways motions there must naturally be a hori- 
zontal pressure gradient which would also act in the ver- 
tical direction and slow down the ascent. This mechanism 
is known in compressible convection as buoyancy braking 
(Hurlburt et al. 1984). 

In order to clarify the nature of the additional braking 
effect seen in the simulations, we first compare with a 
simulation using a somewhat taller box to see whether or 
not the braking sets in later, as would be the case if the 
impenetrable top boundary was the reason for the braking 
effect. (In the following we use calculations where Rq = 1.) 

3.2. Moving the top boundary further away 

With ^0 = 0.3 and ztop = —6.0 we have Zqo = —6.75, see 
Eq. (H) , so we can move the top boundary upwards by no 
more than about 10%. In the following, we discuss a model 
with Ztop = —6.5, but otherwise the same stratification. 
This means that at the new boundary we have to change 
^0 by = 0.4 X Az = 0.2, so we have to require = 0.1. 



Fig. 7. Height and speed of the bubble shown in Fig. ||. 
F = 0.005, 50^ meshpoints. As ^ 0.5, Ro = 0.5. 




Fig. 8. Momentum and integrated buoyancy force acting 
on the bubble shown in Fig. |[ = 0.005, 50'^ meshpoints. 
As = 0.5, Ro = 0.5. 



cause of pressure build-up between the bubble and the top 
boundary. A second possibility could be that the bubble 
loses momentum due to some genuine resistance mecha- 
nism, such as wave braking or viscous friction. However, 
the braking effect seen in the simulations is too sudden 
and too strong to be explained by any genuine braking 
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Fig. 9. Vertical profile of s along the axis of the bubble for 
t = 11 and 2top — —6.5 and = 0.1 (solid line) and ztop — 
—6.0 and — 0.3 (dashed and dash-dotted lines for closed 
and open boundaries). Note that the entropy profile at 
the location of the bubble is not significantly affected by 
the value of ztop- The region where s > Scrit is shown in 
grey. The entropy drop near the surface is a consequence of 
fixing the top temperature, but the location of this entropy 
drop moves further away as we extend the box. F = 0.001, 
50^ X 100 meshpoints, Ro = 1. 

In Fig. ^ we compare the vertical entropy profiles along 
the axis of the bubble for ztop — —6.5 (with = 0.1) 
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Fig. 10. Grey scale representation of the pressure fluctuation togetlier with velocity vectors and a contour marking 
the location where s = Scrit- Light refers to high pressure fluctuation and dark to low pressure fluctuation. The surface 
at z = —6 is marked by a dash-dotted line. Above this line there is a zero-gravity 'buffer layer', modelling the effects 
of an open boundary. 



and Ztop = —6.0 (with = 0.3). We also compare with 
the case of an open top boundary that we have modelled 
by putting an extra layer on top of the box where gravity 
goes smoothly to zero and, as in Brandenburg et al. (1996), 
radiative diffusion is replaced by a heating/cooling term of 
the form —T^^(z)p(T — Ttop), where = everywhere 
except above z — ztop where it goes smoothly to — 
10. This procedure allows the flow to penetrate the layer 
z = Ztop freely. 

It turns out that the entropy profiles at the location 
of the bubble are not significantly affected by the prop- 
erties of the top boundary. The entropy drop near the 
surface is a consequence of fixing the top temperature, 
Ttop. Any increase in the logarithmic pressure at the top, 
(51nptop, causes a corresponding decrease in the entropy, 
Sstop — — 0.4^1nptop- Note, however, that the location of 
this entropy drop at the surface moves further away from 
the location of the bubble as we extend the box. 

In Fig. ^ we show the pressure fluctuations (relative to 
the horizontal mean) together with velocity vectors. Near 
the top of the bubble there is a strong local maximum 
of the pressure fluctuation that drives the gas sideways. 
We have checked that the ram pressure integrated over 
the projected surface is roughly what is needed to explain 
the discrepancy between acceleration and buoyancy force. 
This is suggestive of buoyancy braking being the cause of 
the sudden drop of momentum seen in Fig. ^ (for Rq = 
0.5) and in Fig. |l^ (for the present case of i?o = !)• 

In Fig. [ll] we compare the evolution of h in the two 
cases with different values of Ztop- Within the range of ac- 
curacy the two curves are consistent. Of course, the differ- 
ence in the value of Ztop is not very large, but the increase 
in the total number of scale heights covered in the simula- 



tion is significant: the value of ln(pbot/ptop) has increased 
from 6.2 to 9.0 pressure scale heights. 
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Fig. 11. Height of the bubble (as measured by /i, and 
/itop) for different values of Ztop- Note that the top of the 
bubble, /itop, reaches the top boundary. F = 0.001, 50^ 
and 50^ x 100 meshpoints. As = 0.5, i?o = 1- 

Note that h reaches a maximum around 3, i.e. sig- 
nificantly less than the value of Iztopl- This is partly a 
geometrical effect, because smaller bubbles are able to 
travel somewhat further (in the previous subsection, where 
f?o = 0.5 instead of 1.0, the bubble went to ft. = 3.5). In 
the case of the open top boundary, the bubble rises till 
h = A, which is still small compared with |ztop| = 6. The 
relatively small values of max(/i) are partly due to the fact 
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Fig. 12. The effect of the vertical extent of a weakly 
stratified box on the balance between momentum and in- 
tegrated buoyancy force acting on the bubble. F = 0.005, 
50'^ and 50^ x 100 meshpoints, respectively. As = 0.5, 
Ro = 1. 

that h measures the position of the centre of mass, but in 
the strongly stratified case most of the mass is at the bot- 
tom of the bubble. Furthermore, in all cases considered 
the bubbles take mushroom form with a significant por- 
tion of the mass residing in the stem of the mushroom. A 
somewhat better representation of where most of the hot 
material resides is gained by looking at the value of the 
entropy weighted height, — —z^, where 

= y z{s - s„\t)pdV I J {s ~ Scrit) pdV, (24) 

which is also plotted in Fig. ^ However, even the entropy 
weighted height of the bubble is not very close to the top 
of the bubble. Nevertheless, the location htop of the top 
of the bubble, i.e. where s = Scrit, does reach value close 
to l^topl; see Fig. |ll|. There remains however some worry 
that the centre of mass of the bubble is generally unable 
to travel great distances. In order to clarify this possibility 
we now consider the case of weak stratification where we 
can easily increase the extent of the box. 

3.3. Weak stratification 

We now consider a case of weak stratification and two 
different values of ztop (—6 and —14); hence we choose 
^0 = 30 and — 26.8 respectively, so that the back- 
ground stratification is the same for —6 < z < 2. The 
total number of pressure scale heights in these two cases 
is ln(pbot/ptop) = 0.25 and 0.5, respectively. 

In Fig. |12| we compare the momentum balance for the 
two cases with different values of ztop. Again, the agree- 
ment between the momentum and the time integrated 
buoyancy force is good up to t = 10 (for ztop ~ —6) 
or up to i = 30 (for ztop = —14). The presence of the 



boundary clearly influences the motion of the bubble; nev- 
ertheless the top of the bubble does manage to reach the 
boundary, as we see from Fig. ^ (for weak initial per- 
turbations) and Fig. |l^ (for strong initial perturbations). 
This is in contrast to the results of Sinjab et al. (1990) 
who conclude that the blobs in their calculations 'fail to 
reach the top of the envelope'. In our calculation the hot 
part of the bubble seems to stop at some point just be- 
low the boundary and stays there. Furthermore, Sinjab 
et al. (1990) interpret their two-dimensional calculations 
as indicating a fragmentation. We found no indication of 
any such effect, irrespective of whether the stratification 
was weak or strong. For comparison we also carry out 
two-dimensional cartesian calculations verifying that in 
the two-dimensional case the eddies do eventually travel 
downwards when they come sufficiently close to the hori- 
zontal boundaries. Again, in the weakly stratified case the 
bubbles can travel to large heights, but the effects of the 
boundaries begin to become important much earlier than 
in the three-dimensional case; see Fig. [l^ . 

The results of Sinjab et al. (1990) are of course for the 
stratified case. However, our point is that it is the restric- 
tion to two dimensions (in cartesian geometry) that causes 
boundary effects to become extremely pronounced. Since 
stratification itself can also act as an effective boundary, 
it was necessary to go to the weakly stratified case where 
it is possible to move the boundaries much further away. 
Although the box shown in Fig. |l^ was big enough to 
prevent the edge of the bubble from moving down again, 
boundary effects did begin to affect the evolution at times 
as early as t = 20. 

The ability of our bubble to rise right to the top of 
the box might seem surprising when one recalls that at- 
mospheric thermals only manage to reach a certain maxi- 
mum height. However one must not forget that the buoy- 
ancy force, relative to gravity, is here orders of magnitude 
greater than in the case of a typical atmospheric thermal 
(see e.g. Scorer 1957). A better analogy might therefore be 
that of an air bubble rising through water and eventually 
reaching the surface. 

4. Calculation of the entrainment parameter 

The introduction of a parameter to describe the entrain- 
ment of fiuid by a rising 'cloud' was proposed by Morton 
et al. (1956). They introduced an equation of the form 

V = AirR^a^v, (25) 

which they regarded as describing 'conservation of vol- 
ume'. We have here rewritten their Eq. (16) using a slight 
change of notation. The 'constant' will be referred to 
hereafter as the volume entrainment coefficient. The above 
equation can be simplified to 

R = avW, (26) 
where R is the volume radius. 




Fig. 13. Velocity vectors superimposed on a grey scale representation of the entropy for a tall box with weak strati- 
fication. The velocity is shown in a fixed frame of reference. (The velocity vectors at t = at z « —4 result from the 
initial perturbation.) The initial specific entropy excess of the bubble is As = 0.5 and its initial radius is Ro — 1.0. 
50^ X 100 meshpoints. 



t= D.O f=10.Q i=£0.0 t=40.0 i=7Q.O ^=100.0 




Fig. 14. As for the run shown in Fig. but with strong three-dimensional perturbations. For each time the grey 
scale is here adjusted between minimum and maximum values of s. Note that the bubble still makes it all the way to 
the top of the box, albeit at a somewhat later time. 



We note that Eq. (|2^) is taken over in the work of 
Turner (1963), with becoming Turner's 'alpha'j^ 

Now the concept of 'conservation of volume' lying be- 
hind Eq. (|2^) will be unfamiliar to many physicists. We 
therefore thought it to be worthwhile to concentrate in- 



^ We note that although a has the same meaning in Morton 
et al. (1956) and Turner (1963) the characteristic velocities of 
the bubbles are defined differently, and this is 'compensated' 
by the inclusion of a fc-factor in the entrainment equation 
of Morton et al. (fc — ratio of mean to axial velocity). This 
would not matter, except that when comparing with experi- 
ment, Morton et al. believe the observations relate to the axial 
velocity whereas Turner takes them as referring to the mean 
velocity. This means that when comparing experimental results 
it is really the ak of Morton et al. which should be compared 
with the a of Turner. 



stead on the accretion of mass rather than of volume and 
to write 

M ^ ATrR'^{p)amV, (27) 

where M is the mass entrainment rate, (p) some average 
of the surrounding material density near the bubble, and 
ttni the mass entrainment coefficient. Although an average 
of p over the bubble surface might be more appropriate, 
it is more convenient (and adequate for present purposes) 
to use a volume average. This gives 

M = ^Mda^nV (28) 

where, as before, is the mass displaced by the bubble. 

We can now use Eqs (|6|) and (|2|) to calculate and 
am, all other quantities being already known. The results 
are shown in Figs |l^ and ^ 




Fig. 15. Velocity vectors superimposed on a grey scale representation of the entropy for the same case as in Fig. g_3|, 
but for a two-dimensional cartesian calculation. After t = 20 — 30 the results begin to be affected by boundary effects. 
This is related to the strong nonlocality of two-dimensional calculations. Note that Lx — 16, but in the first two panels 
only a smaller range is shown. i?o = 1.0. 200 x 100 meshpoints. 




5 10 15 20 25 30 indication of a plateau in a^a and a^. Rq = 1.0. 

t 



Fig. 16. Entrainment parameters am (solid line) and 
(dashed line) for the case of strong stratification, ztop = 
—6.0 and — 0.3, and two different initial bubble radii. 
In the upper panel, thick lines (solid and dashed) refer to 
the case with ztop = —6.5 and = 0.1. Note that in the 
case of open boundaries (and ztop = —6.0; upper panel) 
the value of am remains at about 0.4 after t — 10. 



It is interesting to note that for those parts of the 
curves not influenced by the special effects discussed pre- 
viously, the quantity is in reasonable agreement with 
the experimental values of 0.25-0.34 (Morton et al. 1956, 
Turner 1963), with the agreement being better in the 



more strongly stratified case. For smaller initial bub- 
ble radii, Rq — 0.5, attains noticeably larger values. 
Nevertheless, for both bubble radii investigated am is 
roughly comparable (am ~ 0.25). Thus, the assumption 
of Turner (1963) that the value of Rq does not influence 
'alpha' can be verified only for am. 

Better agreement between the 'experimental' and 'cal- 
culated' values cannot perhaps be expected since the 
Reynolds numbers for our calculations are lower than 
those appropriate to the experimental results, which leads 
to a more laminar entrainment process in our case. 
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Regarding am, we note that the entrainment of mate- 
rial leads to an effective drag force on the bubble given 

by 



Drag force = Mv ~ Att pa^^v'^ ^ 



(29) 



which may be compared with the normal hydrodynamic 
drag 



Hydrodynamic drag = ^nR^ pCov 



(30) 



with Cd being the hydrodynamic drag coefficient. 
According to Moore (1967), normal hydrodynamic drag 
plays only a minor role in influencing the motion of a 
thermal. We can now check Moore's point with the help 
of Eqs. 



( |29| ) and (30); we confirm that hydrodynamic drag 
really does have only a minor influence on the motion, ex- 
cept possibly during the final stages in the weakly strati- 
fied case. 

5. Details of the entrainment 

Details of the entrainment process are viewed best in 
terms of tracer particles that are passively advected by the 
flow. In Fig. |l^ we show the location of initially uniformly 
distributed particles at t — 10. Particles that were orig- 
inally inside the bubble (as defined by s > Scrit) remain 
inside the bubble for all times. However, an increasing por- 
tion of new particles from outside the bubble is constantly 
being entrained, which happens mostly through the top 
boundary of the bubble. These particles then move along 
the periphery of the bubble tailwards where they find their 
way into the centre of the ring vortex associated with the 
bubble. 

In Fig. |l^ we have also indicated the trajectory of three 
neighbouring particles that were originally above the bub- 
ble and were subsequently entrained and lifted upwards 
together with the bubble. The kinks in each of these tra- 
jectories correspond to the moment at which the particles 
were overtaken by the bubble and became entrained. 

6. Conclusions 

The main aim of this paper was to test the validity of the 
bubble or thermal concept; in this respect the following 
conclusions can be drawn. 

We found that the bubble could easily be followed 
as a well-defined entity throughout the calculations. 
Nevertheless, its dynamical behaviour consisted of two dis- 
tinct phases. In the first of these, the dynamical behaviour 
expected from the simple bubble dynamics of equating 
buoyancy forces to rate of momentum change was indeed 
confirmed. However, in the second phase an unexpected 
braking effect appeared. Further investigation led us to 
attribute this effect to a combination of boundary effects 
(artificial) and buoyancy braking (real). 

We found, in contrast to the results of Sinjab et al. 
(1990), that the top part of the bubble goes on rising until 
it reaches the surface. Bubbles penetrating to the surface 
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Fig. 18. Tracer particles superimposed on a grey scale 
representation of entropy (bright indicates high entropy). 
Large dots represent particles for which initially s > Scrit , 
i.e. which originated in the initial bubble. Small dots 
represent particles that come from outside the original 
bubble. The contours show the position where s = Scrit 
(dashed for the initial time and solid for t — 10). The 
three neighbouring lines show the particle trajectory of 
entrained particles that were originally outside the bub- 
ble (at X — —0.72, to —0.40 and z = —1.43, and have now 
moved upwards to the point indicated by a small dot). 
Note that — 4, but only the range \x\ < 2.5 is shown. 
100^ X 50 meshpoints. As = 0.5, Rq = 1.0. 

layers favour the view that dissipation is an important 
factor in understanding contact binaries (Hazlehurst 1985, 
1996). 

Since the code used by us was a three-dimensional one 
including effects of viscosity and thermal (radiative) con- 
ductivity, we were in a position to follow in detail and 
more realistically the mass changes of the bubble due to 
entrainment of material. This made it possible for us to 
calculate the value of the entrainment coefficient a rather 
than, as was previously necessary, to take some value from 
glider observations (atmospheric thermals) or experiment. 
The values found by us were in good agreement with the 
those derived empirically, although the non-constancy of 
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'alpha' means that we have throughout referred to the 
entrainment coefBcicnt rather than the entrainment con- 
stant. FinaUy, the introduction of a mass entrainment co- 
efficient to replace (or at least supplement) the previously 
used volume coefficient appears to us desirable on physical 
grounds. 
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